%% 
clc
clear

%% Set initial value
high=4.05;	% height
width=3.3;	% width
Length=20.7;	% Total length
distance=width;	% Wheel spacing

%% Calculate the windward surface area
% Five vehicle model parameters
a1_width=1/2*width;	% TP1
b1_high=1/3*high;
c11=1.5*distance;
c12=distance;
c13=Length-c11-c12;

a2_width=1/2*width;	% TP2
b2_high=high;
c21=2*distance;
c22=0;
c23=Length-c21-c22;

a3_width=1/2*width;	% TP3
b3_high=2/3*high;
c31=2*distance;
c32=distance;
c33=Length-c31-c32;

a4_width=1/2*width;	% TP4
b4_high=1/3*high;
c41=3*distance;
c42=distance;
c43=Length-c41-c42;

a5_width=1/2*width;	% TP5
b5_high=high;
c51=2*distance;
c52=distance; 
c53=Length-c51-c52;


% Part 1 Elliptical Paraboloid
a_width = [a1_width,a2_width,a3_width,a4_width,a5_width]; % TP1-TP5 High speed rail width
b_high = [b1_high,b2_high,b3_high,b4_high,b5_high];  % TP1-TP5 High speed rail height
c = [c11,c21,c31,c41,c51];  % TP1-TP5 High speed rail depth

for i = 1:5   % For loop calculation of TP1-TP5 high-speed railway elliptical paraboloid area
    res(i) = ellipse(a_width(i),b_high(i),c(i));	% Calling the ellipse function to calculate the area res of an elliptical paraboloid
end

s11 = res(1);
s21 = res(2);
s31 = res(3);
s41 = res(4);
s51 = res(5);


% Part 2: Area on the side of the circular frustum
r11=0.5*2/3*high;
s12=pi*(r11+0.5*high)*((0.5*high-r11)^2+c12^2)^0.5;

r21=0;
s22=0;

r31=2/3*high;
s32=0.5*pi*(r31+high)*((high-r31)^2+c32^2)^0.5;

r41=2/3*high;
s42=pi*(r41+high)*((high-r41)^2+c42^2)^0.5;

r51=high;
s52=0.5*pi*high*(Length^2+high^2)^0.5;


% Part III Cylinder Side Area
s13 = 2*pi*0.5*high*c13-width*c13;

s23 = 2*pi*0.5*high*c23-width*c23;

s33 = 2*pi*0.5*high*c33-width*c33;

s43 = 2*pi*0.5*high*c43-width*c43;

s53 = 2*pi*0.5*high*c53-width*c53;


% Total obstructed area of TP1-TP5 high-speed rail
S(1) = s11+s12+s13;
S(2) = s21+s22+s23;
S(3) = s31+s32+s33;
S(4) = s41+s42+s43;
S(5) = s51+s52+s53;


%% Derivative part
syms x y

f = @(x,y)2*x+y*cos(1/3*pi);
diff(f,x,1);	%Taking the first derivative of x
diff(f,y,1);	%Taking the First Derivation of y
F = diff(f,x,y);	%Total differentiation

g = @(x,y)x+y;
diff(g,x,1);	
diff(g,y,1);	
G = diff(g,x,y);	


%% Air resistance iteration (front resistance+air viscosity)

% v_wind=0-13.622       % General situation of wind speed (0-6 levels)��0-49km/h
% v_wind=13.9-28.356     % Special wind speed conditions (7-10 levels)��50-102km/h
% v_rail=0-300       % High-speed rail speed per hour��0-300km/h
% F=1/2*c*p*S*v^2;  % A functional formula for the resistance of the front end to wind

%%====Special weather conditions================
% P=101.3; %101.3KPA
% u=0.5;  %The relative humidity of the air involves 0.8 on rainy days and 0.5 on snowy days;
% t=-10;%10�棬0�棬-10��
% if t==10
%     Pb=1127.9;
% elseif t==0
%     Pb=610.6;
% else
%     Pb=287.98;
% end
% T=273+t;
% p_w=3.48*P/T*(1-0.378*u*0.001*Pb/P);
% p=p_w;
%==============================

afa = 1/3*pi;
p = 1.29;
c = 0.48;
Z = zeros(5,180);	% 5*180
Lj = 25;	% Set one carriage as 25m

t_yiban = 20;	% Under normal circumstances, the temperature is 20 ��
T = t_yiban+273.15;

miu = 1.7894*exp(-5)*(T/288.15)^1.5*((288.15+110.4)/(T+110.4));	% Viscosity level at 15 ��
S_cs = width*Lj*8+high*Lj*8*2;

n=0;
v_rail = 0; %High-speed rail speed per hour
for i = 1:5	%The for loop stores the resistance values of TP1-TP5 high-speed rail into the matrix
    for t = 1:10:1800    %The observation time is 1800s, and the observation values are recorded every 10s
        n = n+1;
        v_wind=0.00074*t;    
        %%v_wind=0.00203*t;   %Overall increase in wind speed under extreme weather conditions
        
        if t>0 && t<=300    %The speed of high-speed rail during different acceleration periods
            v_rail = 0.278*t;
            d1=0.55637;
            d2=0.27874;
        elseif t>300 && t<=1500
            v_rail = 83.4;
            d1=0.00037;
            d2=0.00074;
        else
            v_rail = -0.278*t+500.4;
            d1=-0.55563;
            d2=-0.27726;
        end
        
        V_ct = 2*v_rail+v_wind*cos(afa);	% The relative velocity of air and objects in the front of the vehicle
        V_cs = v_wind+v_rail;	% Relative speed on the body surface
        
        HSR_ct = 0.5*c*p*S(i)*V_ct^2;	% Resistance of high-speed railway locomotive head
        HSR_ctnx = S(i)*miu*d1;
        HSR_csnx = S_cs*miu*d2;
        HSR = HSR_ct+HSR_ctnx+HSR_csnx;
        Z(i,n) = HSR;
        
    end
    n=0;
end
